#include <omp.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include "common.h"

/*
 *  OpenMP version
 */
float *calculation_omp(float *inA, float *inB, int row, int col) {
   float *out;
   int sizeIn, sizeOut;

   if ((row % TILESIZE) != 0) {
      fprintf(stderr,"Number of rows is not a multiple of TILESIZE\n");
      exit(1);
   }
 
   sizeIn = row*col*sizeof(float);
   sizeOut = row*row*sizeof(float);
   out = (float *) malloc(sizeOut);

   int i,j,k;
   
   #pragma omp parallel
   {
      if (omp_get_thread_num() == 0)  
         printf("OpenMP running on %d threads\n",omp_get_num_threads());
   }

   #pragma omp parallel for private(j,k) 
   for (i=0; i < row; i++) {

      float a_x,b_x;
      float sum_a, sum_b, sum_a2, sum_b2, sum_ab, corrcoef;

//     #pragma omp critical
//     fprintf(stderr,"%d,%d\n",omp_get_thread_num(),i);

      for (j=0; j < row; j++) {
         sum_a = sum_a2 = sum_b = sum_b2 = sum_ab = 0;
         for (k=0; k < col; k++) {
            a_x = inA[i*col+k];
            b_x = inB[j*col+k]; 
            sum_a += a_x;
            sum_a2 += a_x*a_x;
            sum_b += b_x;
            sum_b2 += b_x*b_x;
            sum_ab += a_x*b_x;
         }
         corrcoef = (col*sum_ab - sum_a*sum_b)/sqrtf((col*sum_a2-sum_a*sum_a)*(col*sum_b2-sum_b*sum_b));
         out[i*row+j] = corrcoef;
      }

   }
 
   return out;  
}

